Take-home-project

Geospatial Proj

Published

November 27, 2023

Modified

December 2, 2023

Getting Started

Lets make make sure that sfdep, sf, tmap and tidyverse packages of R are currently installed 

pacman::p_load(sf, sfdep, tmap, tidyverse, knitr, plotly, zoo, Kendall, spdep)

Importing OD data into R

Firstly we will import the Passenger volume by Origin Destination Bus Stops data set downloaded from LTA DataMall by using read_csv() of readr package.

odbus <- read_csv("data/aspatial/origin_destination_bus_202308.csv")
glimpse(odbus)
Rows: 5,709,512
Columns: 7
$ YEAR_MONTH          <chr> "2023-08", "2023-08", "2023-08", "2023-08", "2023-…
$ DAY_TYPE            <chr> "WEEKDAY", "WEEKENDS/HOLIDAY", "WEEKENDS/HOLIDAY",…
$ TIME_PER_HOUR       <dbl> 16, 16, 14, 14, 17, 17, 17, 17, 7, 17, 14, 10, 10,…
$ PT_TYPE             <chr> "BUS", "BUS", "BUS", "BUS", "BUS", "BUS", "BUS", "…
$ ORIGIN_PT_CODE      <chr> "04168", "04168", "80119", "80119", "44069", "4406…
$ DESTINATION_PT_CODE <chr> "10051", "10051", "90079", "90079", "17229", "1722…
$ TOTAL_TRIPS         <dbl> 7, 2, 3, 10, 5, 4, 3, 22, 3, 3, 7, 1, 3, 1, 3, 1, …

Origin & Destination Bus Stop Code

odbus$ORIGIN_PT_CODE <-
as.factor(odbus$ORIGIN_PT_CODE)
odbus$DESTINATION_PT_CODE <-
as.factor(odbus$DESTINATION_PT_CODE)

Extracting study data

We will filter out the data according to the requirements set out by Professor 1.) “Weekday @ 6-9am”

origtrip_6_9 <- odbus %>%
  filter(DAY_TYPE == "WEEKDAY") %>%
  filter(TIME_PER_HOUR >= 6 &
           TIME_PER_HOUR <= 9) %>%
  group_by(ORIGIN_PT_CODE) %>%
  summarise(TRIPS = sum(TOTAL_TRIPS))
kable(head(origtrip_6_9))
ORIGIN_PT_CODE TRIPS
01012 1973
01013 952
01019 1789
01029 2561
01039 2938
01059 1651 
write_rds(origtrip_6_9, "data/rds/origtrip_6_9.rds")

2.) “Weekday @ 5-8pm”

origtrip_17_20 <- odbus %>%
  filter(DAY_TYPE == "WEEKDAY") %>%
  filter(TIME_PER_HOUR >= 17 &
           TIME_PER_HOUR <= 20) %>%
  group_by(ORIGIN_PT_CODE) %>%
  summarise(TRIPS = sum(TOTAL_TRIPS))
glimpse(origtrip_17_20)
Rows: 5,035
Columns: 2
$ ORIGIN_PT_CODE <fct> 01012, 01013, 01019, 01029, 01039, 01059, 01109, 01112,…
$ TRIPS          <dbl> 8448, 7328, 3608, 9317, 12937, 2133, 322, 45010, 27233,…
kable(head(origtrip_17_20))
ORIGIN_PT_CODE TRIPS
01012 8448
01013 7328
01019 3608
01029 9317
01039 12937
01059 2133 
write_rds(origtrip_17_20, "data/rds/origtrip_17_20.rds")

3.) “Weekend @ 11am-2pm”

origtrip_11_14 <- odbus %>%
  filter(DAY_TYPE == "WEEKENDS/HOLIDAY") %>%
  filter(TIME_PER_HOUR >= 11 &
           TIME_PER_HOUR <= 14) %>%
  group_by(ORIGIN_PT_CODE) %>%
  summarise(TRIPS = sum(TOTAL_TRIPS))
kable(head(origtrip_11_14))
ORIGIN_PT_CODE TRIPS
01012 2273
01013 1697
01019 1511
01029 3272
01039 5424
01059 1062 
write_rds(origtrip_11_14, "data/rds/origtrip_11_14.rds")

4.) “Weekend @ 4-7pm”

origtrip_16_19 <- odbus %>%
  filter(DAY_TYPE == "WEEKENDS/HOLIDAY") %>%
  filter(TIME_PER_HOUR >= 16 &
           TIME_PER_HOUR <= 19) %>%
  group_by(ORIGIN_PT_CODE) %>%
  summarise(TRIPS = sum(TOTAL_TRIPS))
kable(head(origtrip_16_19))
ORIGIN_PT_CODE TRIPS
01012 3208
01013 2796
01019 1623
01029 4244
01039 7403
01059 1190 
write_rds(origtrip_16_19, "data/rds/origtrip_16_19.rds")
origtrip_11_14 <- read_rds("data/rds/origtrip_11_14.rds")
origtrip_16_19 <- read_rds("data/rds/origtrip_16_19.rds")
origtrip_17_20 <- read_rds("data/rds/origtrip_17_20.rds")
origtrip_6_9 <- read_rds("data/rds/origtrip_6_9.rds")

Importing geospatial data

busstop <- st_read(dsn = "data/geospatial",
                   layer = "BusStop") %>%
  st_transform(crs = 3414)
Reading layer `BusStop' from data source 
  `/Users/youting/ytquek/ISSS624/Project/data/geospatial' using driver `ESRI Shapefile'
Simple feature collection with 5161 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 3970.122 ymin: 26482.1 xmax: 48284.56 ymax: 52983.82
Projected CRS: SVY21

Glimpse the Bus Stop tibble data frame

glimpse(busstop)
Rows: 5,161
Columns: 4
$ BUS_STOP_N <chr> "22069", "32071", "44331", "96081", "11561", "66191", "2338…
$ BUS_ROOF_N <chr> "B06", "B23", "B01", "B05", "B05", "B03", "B02A", "B02", "B…
$ LOC_DESC   <chr> "OPP CEVA LOGISTICS", "AFT TRACK 13", "BLK 239", "GRACE IND…
$ geometry   <POINT [m]> POINT (13576.31 32883.65), POINT (13228.59 44206.38),…

Load Map into MPSZ

mpsz <- st_read(dsn = "data/geospatial",
                   layer = "MPSZ-2019") %>%
  st_transform(crs = 3414)
Reading layer `MPSZ-2019' from data source 
  `/Users/youting/ytquek/ISSS624/Project/data/geospatial' using driver `ESRI Shapefile'
Simple feature collection with 332 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS:  WGS 84
glimpse(mpsz)
Rows: 332
Columns: 7
$ SUBZONE_N  <chr> "MARINA EAST", "INSTITUTION HILL", "ROBERTSON QUAY", "JURON…
$ SUBZONE_C  <chr> "MESZ01", "RVSZ05", "SRSZ01", "WISZ01", "MUSZ02", "MPSZ05",…
$ PLN_AREA_N <chr> "MARINA EAST", "RIVER VALLEY", "SINGAPORE RIVER", "WESTERN …
$ PLN_AREA_C <chr> "ME", "RV", "SR", "WI", "MU", "MP", "WI", "WI", "SI", "SI",…
$ REGION_N   <chr> "CENTRAL REGION", "CENTRAL REGION", "CENTRAL REGION", "WEST…
$ REGION_C   <chr> "CR", "CR", "CR", "WR", "CR", "CR", "WR", "WR", "CR", "CR",…
$ geometry   <MULTIPOLYGON [m]> MULTIPOLYGON (((33222.98 29..., MULTIPOLYGON (…

Geospatial Data Wrangling

Combining Busstop and mpsz

busstop_mpsz <- st_intersection(busstop, mpsz) %>%
  select(BUS_STOP_N, SUBZONE_C) %>%
  st_drop_geometry()
glimpse(busstop_mpsz)
Rows: 5,156
Columns: 2
$ BUS_STOP_N <chr> "13099", "13089", "06151", "13211", "13139", "13109", "1311…
$ SUBZONE_C  <chr> "RVSZ05", "RVSZ05", "SRSZ01", "SRSZ01", "SRSZ01", "SRSZ01",…
write_rds(busstop_mpsz, "data/rds/busstop_mpsz.csv")  
origin_6_9 <- left_join(origtrip_6_9 , busstop_mpsz,
            by = c("ORIGIN_PT_CODE" = "BUS_STOP_N")) %>%
  rename(ORIGIN_BS = ORIGIN_PT_CODE,
         ORIGIN_SZ = SUBZONE_C) %>%
  group_by(ORIGIN_SZ) %>%
  summarise(TOT_TRIPS = sum(TRIPS))
glimpse(origin_6_9)
Rows: 312
Columns: 2
$ ORIGIN_SZ <chr> "AMSZ01", "AMSZ02", "AMSZ03", "AMSZ04", "AMSZ05", "AMSZ06", …
$ TOT_TRIPS <dbl> 116600, 226521, 199437, 114549, 70770, 102390, 23231, 28011,…
origin_17_20 <- left_join(origtrip_17_20 , busstop_mpsz,
            by = c("ORIGIN_PT_CODE" = "BUS_STOP_N")) %>%
  rename(ORIGIN_BS = ORIGIN_PT_CODE,
         ORIGIN_SZ = SUBZONE_C) %>%
  group_by(ORIGIN_SZ) %>%
  summarise(TOT_TRIPS = sum(TRIPS))
origin_11_14 <- left_join(origtrip_11_14 , busstop_mpsz,
            by = c("ORIGIN_PT_CODE" = "BUS_STOP_N")) %>%
  rename(ORIGIN_BS = ORIGIN_PT_CODE,
         ORIGIN_SZ = SUBZONE_C) %>%
  group_by(ORIGIN_SZ) %>%
  summarise(TOT_TRIPS = sum(TRIPS))
origin_16_19 <- left_join(origtrip_16_19 , busstop_mpsz,
            by = c("ORIGIN_PT_CODE" = "BUS_STOP_N")) %>%
  rename(ORIGIN_BS = ORIGIN_PT_CODE,
         ORIGIN_SZ = SUBZONE_C) %>%
  group_by(ORIGIN_SZ) %>%
  summarise(TOT_TRIPS = sum(TRIPS))

Setting up the Hexagon Grid

Drawing the Hexagon Grid

Drawing the hexagon grid over the mpsz map

area_honeycomb_grid = st_make_grid(mpsz, c(500), what = "polygons", square = FALSE)

To sf and add grid ID

honeycomb_grid_sf = st_sf(area_honeycomb_grid) %>%
  mutate(grid_id = 1:length(lengths(area_honeycomb_grid)))
busstop_honeycomb <- st_intersection(honeycomb_grid_sf,busstop) %>%
  select(BUS_STOP_N, grid_id) %>%
  st_drop_geometry()
write_rds(busstop_honeycomb, "data/rds/busstop_honeycomb.csv")
duplicate <- busstop_honeycomb %>%
  group_by_all() %>%
  filter(n()>1) %>%
  ungroup()
busstop_honeycomb <- unique(busstop_honeycomb)

remove grid without value of 0 (i.e. no points in side that grid)

busstop_honeycomb <- busstop_honeycomb %>%
  filter(!is.na(grid_id) & grid_id > 0)

Assign every Bus Stop with a Grid ID

origin_6_9 <- left_join(busstop_honeycomb, origtrip_6_9,
            by = c("BUS_STOP_N" = "ORIGIN_PT_CODE")) 

origin_6_9 <- origin_6_9 %>%
  filter(!is.na(TRIPS) & TRIPS > 0)


origin_17_20 <- left_join(busstop_honeycomb, origtrip_17_20,
            by = c("BUS_STOP_N" = "ORIGIN_PT_CODE")) 

origin_17_20 <- origin_17_20 %>%
  filter(!is.na(TRIPS) & TRIPS > 0)


origin_11_14 <- left_join(busstop_honeycomb, origtrip_11_14,
            by = c("BUS_STOP_N" = "ORIGIN_PT_CODE")) 

origin_11_14 <- origin_11_14 %>%
  filter(!is.na(TRIPS) & TRIPS > 0)


origin_16_19 <- left_join(busstop_honeycomb, origtrip_16_19,
            by = c("BUS_STOP_N" = "ORIGIN_PT_CODE")) 

origin_16_19 <- origin_16_19 %>%
  filter(!is.na(TRIPS) & TRIPS > 0)

Choropleth Visualisation

Weekday Morning Peak 6am-9am

total_trips_by_grid_6_9 <- origin_6_9 %>%
  group_by(grid_id) %>%
  summarise(total_trips = sum(TRIPS, na.rm = TRUE))

Merge geospatial data

total_trips_by_grid_6_9 <- total_trips_by_grid_6_9 %>%
  left_join(honeycomb_grid_sf, by = c("grid_id" = "grid_id"))

total_trips_by_grid_sf_6_9 <- st_sf(total_trips_by_grid_6_9)

Plot the Choropleth map

tmap_mode("view")

tm_shape(total_trips_by_grid_sf_6_9) +
  tm_fill(
    col = "total_trips",
    palette = "Reds",
    style = "cont",
    title = "Total Trips Taken - Weekday Morning Peak 6-9am",
    id = "grid_id",
    showNA = FALSE,
    alpha = 0.6,
    popup.vars = c(
      "Number of trips: " = "total_trips"
    ),
    popup.format = list(
      total_trips = list(format = "f", digits = 0)
    )
  ) +
  tm_borders(col = "grey40", lwd = 0.4)

Weekday Afternoon Peak 5pm - 8pm

total_trips_by_grid_17_20 <- origin_17_20 %>%
  group_by(grid_id) %>%
  summarise(total_trips = sum(TRIPS, na.rm = TRUE))

Merge geospatial data

total_trips_by_grid_17_20 <- total_trips_by_grid_17_20 %>%
  left_join(honeycomb_grid_sf, by = c("grid_id" = "grid_id"))

total_trips_by_grid_sf_17_20 <- st_sf(total_trips_by_grid_17_20)

Plot the Choropleth map

tmap_mode("view")

tm_shape(total_trips_by_grid_sf_17_20) +
  tm_fill(
    col = "total_trips",
    palette = "Reds",
    style = "cont",
    title = "Total Trips Taken - Weekday Afternoon Peak 5 - 8 pm",
    id = "grid_id",
    showNA = FALSE,
    alpha = 0.6,
    popup.vars = c(
      "Number of trips: " = "total_trips"
    ),
    popup.format = list(
      total_trips = list(format = "f", digits = 0)
    )
  ) +
  tm_borders(col = "grey40", lwd = 0.4)

Weekday/Weekend Morning Peak 11am-2pm

total_trips_by_grid <- origin_11_14 %>%
  group_by(grid_id) %>%
  summarise(total_trips = sum(TRIPS, na.rm = TRUE))

Merge geospatial data

total_trips_by_grid <- total_trips_by_grid %>%
  left_join(honeycomb_grid_sf, by = c("grid_id" = "grid_id"))

total_trips_by_grid_sf <- st_sf(total_trips_by_grid)

Plot the Choropleth map

tmap_mode("view")

tm_shape(total_trips_by_grid_sf) +
  tm_fill(
    col = "total_trips",
    palette = "Reds",
    style = "cont",
    title = "Total Trips Taken - Weekday Morning Peak 6-9am",
    id = "grid_id",
    showNA = FALSE,
    alpha = 0.6,
    popup.vars = c(
      "Number of trips: " = "total_trips"
    ),
    popup.format = list(
      total_trips = list(format = "f", digits = 0)
    )
  ) +
  tm_borders(col = "grey40", lwd = 0.4)

Weekend/Holiday Peak 4pm-7pm

total_trips_by_grid <- origin_16_19 %>%
  group_by(grid_id) %>%
  summarise(total_trips = sum(TRIPS, na.rm = TRUE))

Merge geospatial data

total_trips_by_grid <- total_trips_by_grid %>%
  left_join(honeycomb_grid_sf, by = c("grid_id" = "grid_id"))

total_trips_by_grid_sf <- st_sf(total_trips_by_grid)

Plot the Choropleth map

tmap_mode("view")

tm_shape(total_trips_by_grid_sf) +
  tm_fill(
    col = "total_trips",
    palette = "Reds",
    style = "cont",
    title = "Total Trips Taken - Weekday Morning Peak 6-9am",
    id = "grid_id",
    showNA = FALSE,
    alpha = 0.6,
    popup.vars = c(
      "Number of trips: " = "total_trips"
    ),
    popup.format = list(
      total_trips = list(format = "f", digits = 0)
    )
  ) +
  tm_borders(col = "grey40", lwd = 0.4)

Local Indicators of Spatial Association (LISA) Analysis

Computing Contiguity Spatial Weights

a z-score, a pseudo p-value, and a code representing the cluster type for each statistically significant feature. The z-scores and pseudo p-values represent the statistical significance of the computed index values.

Step 1 Check for null Neighbours

# Check for empty neighbor sets
empty_neighbors <- which(sapply(total_trips_by_grid_sf, length) == 0)

# Print the indices of observations with no neighbors
if (length(empty_neighbors) > 0) {
  cat("Observations with no neighbors:", empty_neighbors, "\n")
} else {
  cat("All observations have at least one neighbor.\n")
}
All observations have at least one neighbor.

Step 2 Create a list of all neighbours

# Create neighbor list
wm08 <- total_trips_by_grid_sf %>%
  mutate(nb = st_contiguity(area_honeycomb_grid))

# Filter out observations where 'nb' contains the value 0
# Assuming '0' is an unwanted value in the 'nb' list
wm08 <- wm08 %>%
  filter(!sapply(nb, function(x) any(x == 0)))

# Now, you can proceed with creating the weights
wm08 <- wm08 %>%
  mutate(
    wt = st_weights(nb, style = "W"),
    .before = 1
  )

Step 3 Plot a map based on the list of neighbours

# Set map to static
tmap_mode("view")

map_08 <- tm_shape(wm08) +
  tm_fill(
    col = "total_trips",
    palette = "PuRd",
    style = "pretty",
    title = "Trip sizes"
  ) +
  tm_layout(main.title = "Total Bus Trips Across the Island ",
            main.title.size = 1,
            main.title.position = "center",
            legend.position = c("left", "top"),
            legend.height = .6,
            legend.width = .2,
            frame = FALSE
            
  )

map_08

Weekday Morning Peak Traffic

# Set seed to ensure that results from simulation are reproducible
set.seed(1234)
# Step 1: Merge the data
# Ensure that both wm08 and Origin_6_9 have 'grid_id' as a common key
# and that it's of the same data type in both data frames
merged_data <- wm08 %>% 
  left_join(origin_6_9, by = "grid_id")

# Step 2: Prepare the data
# Assuming 'area_honeycomb_grid' is the geometry column in wm08
# Convert to an sf object if not already
if (!("sf" %in% class(merged_data))) {
  merged_data <- st_as_sf(merged_data, geom_col = "area_honeycomb_grid")
}

# Recreate the neighbor list and spatial weights
# Ensure 'nb' is in the correct format (e.g., created using st_contiguity or similar)
listw <- nb2listw(merged_data$nb, style = "W")

merged_data$standardized_trips <- as.numeric(scale(merged_data$TRIPS, center = TRUE, scale = TRUE))

# Remove NA values from the data
merged_data <- merged_data[!is.na(merged_data$standardized_trips), ]

# Recreate the spatial weights to match the filtered data
listw <- nb2listw(merged_data$nb, style = "W")

# Check if the lengths match
if (nrow(merged_data) != length(listw$neighbours)) {
  stop("The length of the data and the spatial weights list do not match.")
}

# Run the Monte Carlo test for Local Moran's 
lisa_wd_morn <- merged_data %>%
  mutate(local_moran = local_moran(
    total_trips, nb, wt, nsim = 99),
    .before = 1) %>%
  unnest(local_moran)

# unnest the dataframe column
tidyr::unnest(lisa_wd_morn)
Simple feature collection with 23878 features and 19 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 4167.538 ymin: 27151.39 xmax: 46917.54 ymax: 50245.4
Projected CRS: SVY21 / Singapore TM
# A tibble: 23,878 × 20
      ii      eii var_ii  z_ii  p_ii p_ii_sim p_folded_sim skewness kurtosis
   <dbl>    <dbl>  <dbl> <dbl> <dbl>    <dbl>        <dbl>    <dbl>    <dbl>
 1 0.285 -0.00749 0.0694 1.11  0.267     0.02         0.01    -2.06     5.30
 2 0.285 -0.00749 0.0694 1.11  0.267     0.02         0.01    -2.06     5.30
 3 0.285 -0.00749 0.0694 1.11  0.267     0.02         0.01    -2.06     5.30
 4 0.284  0.0734  0.0609 0.852 0.394     0.02         0.01    -2.12     4.77
 5 0.284  0.00560 0.312  0.498 0.619     0.06         0.03    -4.62    23.5 
 6 0.270  0.00300 0.0768 0.962 0.336     0.04         0.02    -3.85    22.0 
 7 0.270  0.00300 0.0768 0.962 0.336     0.04         0.02    -3.85    22.0 
 8 0.270  0.00300 0.0768 0.962 0.336     0.04         0.02    -3.85    22.0 
 9 0.278  0.0170  0.0599 1.07  0.287     0.02         0.01    -2.71     9.88
10 0.278  0.0170  0.0599 1.07  0.287     0.02         0.01    -2.71     9.88
# ℹ 23,868 more rows
# ℹ 11 more variables: mean <fct>, median <fct>, pysal <fct>, wt <dbl>,
#   grid_id <int>, total_trips <dbl>, area_honeycomb_grid <POLYGON [m]>,
#   nb <int>, BUS_STOP_N <chr>, TRIPS <dbl>, standardized_trips <dbl>

Weekday Afternoon/Evening Peak Traffic

# Set seed to ensure that results from simulation are reproducible
set.seed(1234)
# Step 1: Merge the data
# Ensure that both wm08 and Origin_6_9 have 'grid_id' as a common key
# and that it's of the same data type in both data frames
merged_data_17_20 <- wm08 %>% 
  left_join(origin_17_20, by = "grid_id")

# Step 2: Prepare the data
# Assuming 'area_honeycomb_grid' is the geometry column in wm08
# Convert to an sf object if not already
if (!("sf" %in% class(merged_data_17_20))) {
  merged_data_17_20 <- st_as_sf(merged_data_17_20, geom_col = "area_honeycomb_grid")
}

# Recreate the neighbor list and spatial weights
# Ensure 'nb' is in the correct format (e.g., created using st_contiguity or similar)
listw <- nb2listw(merged_data_17_20$nb, style = "W")

merged_data_17_20$standardized_trips <- as.numeric(scale(merged_data_17_20$TRIPS, center = TRUE, scale = TRUE))

# Remove NA values from the data
merged_data_17_20 <- merged_data_17_20[!is.na(merged_data_17_20$standardized_trips), ]

# Recreate the spatial weights to match the filtered data
listw <- nb2listw(merged_data_17_20$nb, style = "W")

# Check if the lengths match
if (nrow(merged_data_17_20) != length(listw$neighbours)) {
  stop("The length of the data and the spatial weights list do not match.")
}

# Run the Monte Carlo test for Local Moran's 
lisa <- merged_data_17_20 %>%
  mutate(local_moran = local_moran(
    total_trips, nb, wt, nsim = 99),
    .before = 1) %>%
  unnest(local_moran)

# unnest the dataframe column
tidyr::unnest(lisa)
Simple feature collection with 23931 features and 19 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 4167.538 ymin: 27151.39 xmax: 46917.54 ymax: 50245.4
Projected CRS: SVY21 / Singapore TM
# A tibble: 23,931 × 20
      ii      eii var_ii  z_ii  p_ii p_ii_sim p_folded_sim skewness kurtosis
   <dbl>    <dbl>  <dbl> <dbl> <dbl>    <dbl>        <dbl>    <dbl>    <dbl>
 1 0.284  0.0197  0.0964 0.850 0.395     0.02         0.01    -3.57     16.0
 2 0.284  0.0197  0.0964 0.850 0.395     0.02         0.01    -3.57     16.0
 3 0.284  0.0197  0.0964 0.850 0.395     0.02         0.01    -3.57     16.0
 4 0.283  0.0310  0.161  0.628 0.530     0.04         0.02    -4.53     27.4
 5 0.283  0.0399  0.154  0.618 0.537     0.06         0.03    -4.76     30.1
 6 0.269 -0.0326  0.249  0.604 0.546     0.04         0.02    -4.29     20.0
 7 0.269 -0.0326  0.249  0.604 0.546     0.04         0.02    -4.29     20.0
 8 0.269 -0.0326  0.249  0.604 0.546     0.04         0.02    -4.29     20.0
 9 0.277  0.00684 0.0737 0.995 0.320     0.02         0.01    -3.22     13.1
10 0.277  0.00684 0.0737 0.995 0.320     0.02         0.01    -3.22     13.1
# ℹ 23,921 more rows
# ℹ 11 more variables: mean <fct>, median <fct>, pysal <fct>, wt <dbl>,
#   grid_id <int>, total_trips <dbl>, area_honeycomb_grid <POLYGON [m]>,
#   nb <int>, BUS_STOP_N <chr>, TRIPS <dbl>, standardized_trips <dbl>
# Extracting p-values or other results from local_moran_mc
# ...

Weekend/Holiday Morning Peak Traffic

# Set seed to ensure that results from simulation are reproducible
set.seed(1234)
# Step 1: Merge the data
# Ensure that both wm08 and Origin_6_9 have 'grid_id' as a common key
# and that it's of the same data type in both data frames
merged_data_11_14 <- wm08 %>% 
  left_join(origin_11_14, by = "grid_id")

# Step 2: Prepare the data
# Assuming 'area_honeycomb_grid' is the geometry column in wm08
# Convert to an sf object if not already
if (!("sf" %in% class(merged_data_11_14))) {
  merged_data_11_14 <- st_as_sf(merged_data_11_14, geom_col = "area_honeycomb_grid")
}

# Recreate the neighbor list and spatial weights
# Ensure 'nb' is in the correct format (e.g., created using st_contiguity or similar)
listw <- nb2listw(merged_data_11_14$nb, style = "W")

merged_data_11_14$standardized_trips <- as.numeric(scale(merged_data_11_14$TRIPS, center = TRUE, scale = TRUE))

# Remove NA values from the data
merged_data_11_14 <- merged_data_11_14[!is.na(merged_data_11_14$standardized_trips), ]

# Recreate the spatial weights to match the filtered data
listw <- nb2listw(merged_data_11_14$nb, style = "W")

# Check if the lengths match
if (nrow(merged_data_11_14) != length(listw$neighbours)) {
  stop("The length of the data and the spatial weights list do not match.")
}

# Run the Monte Carlo test for Local Moran's 
lisa <- merged_data_11_14 %>%
  mutate(local_moran = local_moran(
    total_trips, nb, wt, nsim = 99),
    .before = 1) %>%
  unnest(local_moran)

# unnest the dataframe column
tidyr::unnest(lisa)
Simple feature collection with 23830 features and 19 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 4167.538 ymin: 27151.39 xmax: 46917.54 ymax: 50245.4
Projected CRS: SVY21 / Singapore TM
# A tibble: 23,830 × 20
      ii     eii var_ii  z_ii  p_ii p_ii_sim p_folded_sim skewness kurtosis
   <dbl>   <dbl>  <dbl> <dbl> <dbl>    <dbl>        <dbl>    <dbl>    <dbl>
 1 0.285 -0.0458 0.154  0.845 0.398     0.02         0.01    -2.90    10.4 
 2 0.285 -0.0458 0.154  0.845 0.398     0.02         0.01    -2.90    10.4 
 3 0.285 -0.0458 0.154  0.845 0.398     0.02         0.01    -2.90    10.4 
 4 0.284  0.0145 0.237  0.554 0.579     0.04         0.02    -4.81    26.3 
 5 0.284  0.0439 0.0833 0.833 0.405     0.02         0.01    -2.43     6.68
 6 0.270  0.0165 0.0608 1.03  0.303     0.04         0.02    -2.18     5.62
 7 0.270  0.0165 0.0608 1.03  0.303     0.04         0.02    -2.18     5.62
 8 0.270  0.0165 0.0608 1.03  0.303     0.04         0.02    -2.18     5.62
 9 0.278 -0.0252 0.0940 0.990 0.322     0.02         0.01    -2.44     6.09
10 0.278 -0.0252 0.0940 0.990 0.322     0.02         0.01    -2.44     6.09
# ℹ 23,820 more rows
# ℹ 11 more variables: mean <fct>, median <fct>, pysal <fct>, wt <dbl>,
#   grid_id <int>, total_trips <dbl>, area_honeycomb_grid <POLYGON [m]>,
#   nb <int>, BUS_STOP_N <chr>, TRIPS <dbl>, standardized_trips <dbl>
# Extracting p-values or other results from local_moran_mc
# ...

Weekend/Holiday Evening Peak Traffic

# Set seed to ensure that results from simulation are reproducible
set.seed(1234)
# Step 1: Merge the data
# Ensure that both wm08 and Origin_6_9 have 'grid_id' as a common key
# and that it's of the same data type in both data frames
merged_data_16_19 <- wm08 %>% 
  left_join(origin_16_19, by = "grid_id")

# Step 2: Prepare the data
# Assuming 'area_honeycomb_grid' is the geometry column in wm08
# Convert to an sf object if not already
if (!("sf" %in% class(merged_data_16_19))) {
  merged_data_16_19 <- st_as_sf(merged_data_16_19, geom_col = "area_honeycomb_grid")
}

# Recreate the neighbor list and spatial weights
# Ensure 'nb' is in the correct format (e.g., created using st_contiguity or similar)
listw <- nb2listw(merged_data_16_19$nb, style = "W")

merged_data_16_19$standardized_trips <- as.numeric(scale(merged_data_16_19$TRIPS, center = TRUE, scale = TRUE))

# Remove NA values from the data
merged_data_16_19 <- merged_data_16_19[!is.na(merged_data_16_19$standardized_trips), ]

# Recreate the spatial weights to match the filtered data
listw <- nb2listw(merged_data_16_19$nb, style = "W")

# Check if the lengths match
if (nrow(merged_data_16_19) != length(listw$neighbours)) {
  stop("The length of the data and the spatial weights list do not match.")
}

# Run the Monte Carlo test for Local Moran's 
lisa <- merged_data_16_19 %>%
  mutate(local_moran = local_moran(
    total_trips, nb, wt, nsim = 99),
    .before = 1) %>%
  unnest(local_moran)

# unnest the dataframe column
tidyr::unnest(lisa)
Simple feature collection with 23843 features and 19 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 4167.538 ymin: 27151.39 xmax: 46917.54 ymax: 50245.4
Projected CRS: SVY21 / Singapore TM
# A tibble: 23,843 × 20
      ii     eii var_ii  z_ii  p_ii p_ii_sim p_folded_sim skewness kurtosis
   <dbl>   <dbl>  <dbl> <dbl> <dbl>    <dbl>        <dbl>    <dbl>    <dbl>
 1 0.285  0.0391 0.0355 1.31  0.192     0.02         0.01    -1.56     2.78
 2 0.285  0.0391 0.0355 1.31  0.192     0.02         0.01    -1.56     2.78
 3 0.285  0.0391 0.0355 1.31  0.192     0.02         0.01    -1.56     2.78
 4 0.284  0.0645 0.0767 0.793 0.428     0.12         0.06    -2.68    10.2 
 5 0.284  0.0481 0.0861 0.804 0.421     0.1          0.05    -3.32    14.7 
 6 0.270 -0.0307 0.0930 0.986 0.324     0.04         0.02    -2.86    11.5 
 7 0.270 -0.0307 0.0930 0.986 0.324     0.04         0.02    -2.86    11.5 
 8 0.270 -0.0307 0.0930 0.986 0.324     0.04         0.02    -2.86    11.5 
 9 0.278  0.0350 0.0613 0.982 0.326     0.02         0.01    -2.95    10.9 
10 0.278  0.0350 0.0613 0.982 0.326     0.02         0.01    -2.95    10.9 
# ℹ 23,833 more rows
# ℹ 11 more variables: mean <fct>, median <fct>, pysal <fct>, wt <dbl>,
#   grid_id <int>, total_trips <dbl>, area_honeycomb_grid <POLYGON [m]>,
#   nb <int>, BUS_STOP_N <chr>, TRIPS <dbl>, standardized_trips <dbl>
# Extracting p-values or other results from local_moran_mc
# ...
# Set map to static
tmap_mode("plot")

# weekday morning moran I values
map_wdm_moran08 <- tm_shape(lisa_wd_morn) +
  tm_fill(
    col = "ii",
    palette = "OrRd",
    style = "pretty",
    title = "Local Moran's I"
  ) +
  tm_layout(main.title = "Weekday Morning Peak Traffic",
            main.title.size = 1,
            main.title.position = "center",
            legend.position = c("left", "top"),
            legend.height = .6,
            legend.width = .2,
            frame = FALSE
  )

map_wdm_moran08

# Set map to static
tmap_mode("view")

map_m1 <- tm_shape(merged_data) +
  tm_fill(
    col = "total_trips",
    palette = "PuRd",
    style = "pretty",
    title = "Trip sizes"
  ) +
  tm_layout(main.title = "Total Bus Trips Across the Island ",
            main.title.size = 1,
            main.title.position = "center",
            legend.position = c("left", "top"),
            legend.height = .6,
            legend.width = .2,
            frame = FALSE
  )

map_m1